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Ab Initio Theory of Superconductivity in a Magnetic Field I. : Spin Density 
Functional Theory For Superconductors and Eliashberg Equations. 
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We present a first-principles approach to describe magnetic and superconducting systems and the 
phenomena of competition between these electronic effects. We develop a density functional theory: 
SpinSCDFT, by extending the Hohenberg-Kohn theorem and constructing the non-interacting Kohn- 
Sham system. An exchange-correlation functional for SpinSCDFT is derived from the Sham Schliiter 
connection between the SpinSCDFT Kohn-Sham and a self-energy in Eliashberg approximation. The 
reference Eliashberg equations for superconductors in the presence of magnetism are also derived 
and discussed. 


I. INTRODUCTION 

In this work, we present how magnetic (M) and super¬ 
conducting (SC) properties can be computed on the same 
footing and from first principles by extending the Density 
Functional Theory (DFT) framework. In developing this 
spin DFT for SC (SpinSCDFT) we will restrict ourselves 
to situations where currents are negligible and only con¬ 
sider the effect of the Zeeman term of the Hamiltonian. 
Under this assumption we can exclude the occurrence of 
the Abrikosov vortex state 1 , that having a mesoscopic 
characteristic length-scale would be beyond the present 
computational power for a fully ab-initio method. 

The expulsion of static M fields from the bulk 2 is one 
of the most spectacular properties of SC materials and 
illustrates the profound competition between M and SC 
behavior. The SC-M interaction generates in fact a large 
number of interesting phenomena on which the scientific 
community has focused attention. Some of the most in¬ 
vestigated are the Abrikosov vortices 1 and the variety of 
fascinating effects occurring in heterostructures 3 , such as 
stacked layers of M with SC material (see Ref. 4 for a re¬ 
view). Among these effects is the FFLO state, named 
after Fulde, Ferrel 5 , Larkin and Ovchinnikov 6 , where 
strong exchange fields induce a SC state with a finite 
momentum pairing. This state was recently observed 
experimentally 7,8 in heavy Fermion SC, many years after 
its prediction. In addition, triplet SC has been observed 
in several systems 9-15 , and is usually associated to ferro¬ 
magnetism. 

Among the many effects generated by the interplay of 
magnetism and superconductivity, some have an intrin¬ 
sic microscopic nature and could be accessible to first- 
principle calculations, in particular we refer to the sharp 
suppression of the critical temperature due to paramag¬ 
netic impurities 16 , and the surprising evidence of coex¬ 
isting phases between singlet SC and local magnetism, 
in particular close to a magnetic phase boundary 14,17-19 
where high— T c SC occurs 20,21 . We devote this work to 
set the ground for an ab-initio theory to describe these 
physical effects. 

We will start our formulation from the Pauli Hamil- 
toninan (Sec. II). In Sec. (Ill), we formulate a density 
functional theory (DFT), proving that the electronic den¬ 


sity n(r) , the spin magnetization m(r), the diagonal of 
the nuclear ./V-body density matrix and the singlet and 
triplet SC order parameters x(r, r') are uniquely con¬ 
nected with their respective external potentials. With 
this extension of the Hohenberg-Kohn theorem 22 we lay 
the foundation of the DFT for M and SC systems: Spin¬ 
SCDFT. In Sec. Ill A we introduce the formally non¬ 
interacting Kohn Sham (KS) system that reproduces the 
exact densities of the interacting system. Similar to ev¬ 
ery DFT, SpinSCDFT relies on the construction of an 
exchange correlation ( xc ) functional that connects the 
KS with the interacting system. In this work, this is 
achieved by establishing, in Sec. IIIB, a Sham-Schliiter 
connection 23 via the Dyson equation of the interacting 
system. 

The interacting system is also being investigated di¬ 
rectly by means of a magnetic extension of the Eliashberg 
method 24-28 . A derivation 29 of this alternative approach 
in the present notation is given in Sec. IV. Advantages 
and disadvantages of these two theoretical schemes, Spin¬ 
SCDFT and Eliashberg, will be discussed in the conclu¬ 
sions. 


II. HAMILTONIAN 

We assume that the interacting system is governed 
by the Pauli Hamiltonian (we use Hartree atomic units 
throughout) 

h = t + r„ + v. + u.e + U en + u nn , (1) 

where T e (T n ) is the kinetic energy operator of the elec¬ 
tron (nuclei) and U ee (U nn ) is the electron-electron (nuclei- 
nuclei) interaction, i.e. usually the Coulomb potential. 
t/ on is the Coulomb potential between electrons and nu¬ 
clei. To break the respective symmetries and allow the 
corresponding densities to adopt non-zero values in a 
thermal average we include an external vector poten¬ 
tial A 0Xt (r) and an external singlet/triplet pair poten¬ 
tial A oxt (r,r') in the Hamiltonian. These external fields 
will be set to zero at the end of the derivation. Because 
we do not consider currents, the only term in the Pauli 
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Hamiltonian containing A oxt (r) is: 

C 2 

f e = jdr ft (r) ■ (-00— + S • B ext (r)) • ^(r) ( 2 ) 

with B ext (r) = V x A oxt (r) and S = §( 0 y o’, ) T , 
<J x ,y z being the Pauli matrices. We use the notation 
ft (r) = ( ft ( r |) ( r |) ) for the field operator where 
ft (r f) creates an electron at location r with spin up. 
The scalar potential part of H reads: 


the KS system is done assuming that densities are always 
v— representable i.e. we assume the existence of the KS 
system. Being non-interacting it consists of independent 
equations for nuclei and electrons, coupled only via the 
xc potentials. Our focus will be on the electronic system, 
discussed in detail in Sec. Ill A 2. The nuclear part will 
be addressed in Sec. IIIA1, briefly, since it is usually 
enough to approximate the nuclear KS system with 
its non SC counterpart 33,35 . The construction the xc 
potentials will be discussed in Sec. IIIB and Sec. IIIC. 


V B = Jdr ft (r) • oo • if{r)v Bxt (r) 

\ f dr f dr '{x( r ’r') ' A““*(r,r') + h.c.) . (3) 
Here, the anomalous density operator is defined by 


X(r, r') = ftr) • $ • ftr') . (4) 


x(r,r') is a 4-vector of which the first component (pro¬ 
portional to $i) is the singlet part of the order param¬ 
eter, while the other components (related to $ 2,$3 and 
4 * 4 ) are the triplet part. The 4 components of the sin¬ 
glet/triplet vector 4? = (io- y , —a z ,ao,a x ) T are 2 x 2 spin 
matrices similar to the components of S. Similarly, the 
anomalous external potential 


/ Ag xt (r, r') 
Z\“‘(r,r') 


( 5 ) 


is assumed to have singlet and triplet components. 


III. SPIN SCDFT 

The conventional density functional approach to the 
Many-Body problem 22, 30-32 consists of two steps: first 
establishing the Hohenberg Kohn (HI<) theorem, i.e. re¬ 
alize that a chosen set of densities is uniquely connected 
with a set of external potentials; second, construct an 
auxiliary, non-interacting KS system to reproduce the 
densities of the interacting system. 

We follow Ref. 33 and consider a multi-component 
DFT with the normal n(r ), the SC order parameter 
as the anomalous density x( r i r '), that describes the 
electrons condensed into singlet and triplet states, and 
r(Ri..R]y ) the diagonal of the nuclear iV-body density 
matrix. In addition, we introduce the magnetization 
m(r) as another electronic density. 

The HK proof (n(r),m(r),x(r,r'),r(Ri..RN)) o 
(iwMj A ext (r, r'), W b ^(Ri..Rn)) is a 

straightforward generalization of Mermin’s HK proof in 
a finite temperature ensemble 34 . For this reason we will 
not repeat it here. On the other hand the construction of 


A. The Kohn-Sham System 

In this work we are mainly interested in the influence 
of a magnetic field on the SC state. We briefly review 
the approximation steps to arrive at the Frohlich Hamil¬ 
tonian starting from the formally exact multi-component 
DFT. The reader may refer to the existing literature for 
further details 32,33 . We introduce the KS Hamiltonian 

H ks = H^s + Hks > ( 6 ) 

where we have separated the electronic 

His = Jdr ft (r) ■ o- 0 (—— + v B {r) -ft) ■ ftr) 

— - jdr Jdr 1 (x(r, r •’) ■ (r, r') + h.c.^ 

+ Jdr rh(r) ■ B s (r), (7) 

from the nuclear 77“ s 

£ks = - JdRC\R)^C(R) + 

+ j-jdRi■ -dR Nn C t ( J Ri).--C t (R Nn ) X 
xW B (fi 1 ,...,2?^JC(i?i)...C(i?ArJ. (8) 

We write v a (r) = u ext (r) + u xc (r) with u xc (r) being the 
scalar xc potential (similar for B s (r) and A B (r,r')). 
rh(r) = ft(r) • S • ftr ) is the operator of the mag¬ 
netic density. In the nuclear description, C (R) cre¬ 
ates the nuclear field at location R. Following Lriders 
et al. 33 and Marques et al. 35 we use the N— body po¬ 
tential W.(Ri,R n J because in this way the nuclear 
KS system can be easily related to the standard Born- 
Oppenheimer approximation. M refers to the ionic mass. 
Here, we neglect the spin of the nuclei and consider only 
one atomic type (the generalization is straightforward). 


1. The Nuclear Part 

Since SC occurs in the solid phase, we assume that ions 
can only perform small oscillations about their equilib¬ 
rium position. A discussion that goes beyond this sim¬ 
ple picture can be found in Ref. 36 and 32. We expand 
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W a (Ri, ..., Rn„) in u.i = Ri — R a j around the equilib¬ 
rium positions Roi- The nuclear degrees of freedom (up 
to harmonic order) are described by the Hamiltonian H £g 
with H™ s = H + 0(u 3 ) in second quantization 

Hf s = J2^{b% + \)- (9) 
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We use the notation q = q, A with Bloch vector q and 
mode number A. We further use the notation —q = q , A 
for all Bloch vector and band or mode combinations. 
We point out that via the functional dependence of 
W s [n, m , Xi T] the KS phonon frequencies fi q are in prin¬ 
ciple functionals of the densities as well. w q creates a 
bosonic KS phonon with quantum numbers q. Usually, 
approximating W B with the Born-Oppenheimer energy 
surface, leads to phonon frequencies in excellent agree¬ 
ment with experiment 37,38 . 

The electron phonon scattering should be formally con¬ 
structed from the bare Coulomb interaction 36 . However 
in order to have a proper description of the electronic 
screening this is not feasible in practice. The solution is 
the substitution of the many body electron phonon in¬ 
teraction with its Kohn Sham counterpart U Bn —► “ ph . 


^Ks ph = f dr 9q(r)^{ r ) • ■ i>(r)(b q + tf_ q ) , ( 10 ) 

q m 


where m = 0,z and g° q {r) = g z q (r) = 5 ^ (r) , u 

being the phononic displacement vectors 3 '’ 38 . This form 
incorporates most of the electronic influence on the bare 
Coulomb interaction between electrons and nuclei. We 
consider this as a good approximation for the dressed 
phonon vertex in the non-SC state, see also Ref. 36 for 
a further discussion. Note that (J?^“ ph ) is part of the 
xc— functional of the electronic KS system and will be 
added later in our approximate functional using pertur¬ 
bation theory. For later use in the derivation of the xc 
potential, we define the propagator of the non-interacting 
system of KS phonons 


KAr) = + h - q to) (MO) + bl q ,(0))U, 

(11) 

( 12 ) 


E q,q'{ V n) ~ &q,-q'^ 


1 


1 


1 ^q 


Here T is the usual time (r) ordering symbol of operators 
& g (r)+ 6 l 9 (r) in the Heisenberg picture and (.. .) ph means 
to evaluate the thermal average using the Hamiltonian 
of Eq. (9). The bosonic Matsubara frequency is 


2. The Electronic Part 

The electronic KS Hamiltonian H^ s is not diagonal in 
the electronic field operator ^(r) because Eq. (7) involves 


terms proportional to ipif) and ifiiffi. Being a hermitian 
operator, we can find an orthonormal set of eigenvectors 
of H ^_ s in which it is diagonal. Let 7 ! create such a two 
component vector in spin space (the Hamiltonian is not 
diagonal in spin so the spin degrees of freedom is in the 
set {fc}), then the SC KS system will take the form 

#k s = E o + ^2 Ek - 0 • (I 3 ) 

k 

where Eq is the ground state energy and the Ek 
are all positive. This form can be achieved 39 by 
commuting the operators H^ s = Ek/j\Ak = 

Yhk\E k <o + Sfc|E fe >o E k(t\b,k + Sfc|B fe <o \ E k\b,ka\ and 
then redefining the negative energy particle operators as 
holes fife = 7 fc- We use a notation that is based on the 
one of Ref. 40, 41 and 24. We introduce 


E(r) 


( \ 

4>(r I) 

p(r f) 

V fo{r 4.) / 


(14) 


Using this Nambu field operator »F(r) the KS Hamilto¬ 
nian reads 

E ls = Jdr Jdr'&(r) ■ ^H KS (r,r') (15) 

where the KS Hamiltonian (first quantization Nambu 
form) is given by 

H KS {r,r') = 

6 (r — r')H£s{r) <& ■ A B (r,r f ) 

-($ • A s (r, r')Y - 6 (r - r')(H™(r)) Ts 

with 

= (^V 2 +v B (r) - M )a 0 - S • B.(r) . (17) 

Note that the changed order of the electronic field oper¬ 
ator implies a transposition in spin space in the (- 1 ,- 1 ) 
component that is equivalent to using S*. In a similar 
transformation the diagonal KS Hamiltonian Eq. (13) be¬ 
comes 

= as) 

k ' ' 



with &k = ^ ^ • As a consequence of the rearrange¬ 

ment of the operators, in the Nambu-Anderson form 
should appear the trace of the Hamiltonian H^ s . How¬ 
ever, not being an operator, this cancels from thermal 
averages and has been disregarded. T>k is a two (not 
four) component vector because the spin may not be a 
good quantum number in the SC KS system. We can 
diagonalize the form in Eq. (15) to the form Eq. (18) 
by introducing a unitary transformation that we param- 
eterise generically with four complex spinor functions. 
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This connection between <T(r) and 9> k is known as the 
Bogoliubov-Valatin transformation 42,43 . We write it in 
the form 


H r ) = Y 

k 


Uk{r) v* k (r) 
Vk(r) u* k ( r ) 


■$k 





u*k( r ) v*k( r ) 
Vk{r) u k (r) 


■ T(r). 


(19) 


Note that in the first case the matrix is 4 x 2 dimensional, 
and in the second 2x4 because of the spinor property of 
the u k (r),v k (r). In going from Eq. (15) to Eq. (18), we 
identify 




(<( r ) ^( r ) 
\Vk{r) u k (r) 


■H K3 {r,r')- 


fuk’ir') v* k ,(r')\ 
V>k'(r') u* k ,lr’)J 


( E k 0 

^ 0 -E k 


&kk' ) 


( 20 ) 


which are the KS Bogoliubov de Gennes (KSBdG) 
equations for magnetic system. Applying the inverse 
Bogoliubov-Valatin transformation from the left we ob¬ 
tain two redundant vector equations of which we usually 
consider the first for the positive eigenvalues E k 

■(£$)-*( $!)■<“) 

This is the usual form of the KSBdG equations which 
generalize those of Ref. 44 and Ref. 33. The equation in 
( ^fc( r ) 'Uk( r ) ) leads to the equivalent negative eigen¬ 
value —E k which reflects the additional degrees of free¬ 
dom that we have created in going to the 2x2 Nambu 
formalism. 

a. The Normal State KS Basis expansion The KS¬ 
BdG equations 21 pose a challenging integro differential 
problem. Sensible approximations can be obtained by 
first performing an expansion into a basis set that is ac¬ 
cessible in practice and resembles closely to the true quasi 
particle structure of the non-superconducting phase of 
the material under consideration. With this in mind we 
consider the non-SC KS single particle equation: 


£i*<Pia(r) = ((-^-+u°(r) -m)oo ~(p ia {r) 

( 22 ) 

v°(r) and B° z (r) are known functionals, like the local 
spin density approximation (LSD A) 45 . We also assume 
that B° s is collinear and has components in cr~ only. We 
use a pure spinor notation for the orbitals, i.e. (pi a (r) 
has only one non-vanishing component, e.g. Ti\{ r ) = 

We use the indices i,j for the quantum 

numbers of the basis and thus distinguish from the quan¬ 
tum number k of the SC KS system. Later, in the Spin 
Decoupling Approximation III A 2 c when we assume the 
expansion coefficients to have only one non-vanishing 
component each, this distinction will not be made. As a 


(* i( o rt) )' 


next step we expand the Bogoliubov-Valatin transforma¬ 
tions in these solutions {</?i CT (r)} 46 

^k( r ) = Y u k , Vk(r) = J2 v k T( fia( r )- (23) 

icr icr 


Defining 

; the matrix elements 



7? crcr ' 

= J dr Tia( cr o(v s (r) - 




-S- {B s (r)-B°J 


(24) 

A eaa ' 

v 

= j dr fdr'&Jr) • ( 

r')) • 

<? 0 Ar') (25) 

ecru' 

c ij 

= Sio-SijSatr' + , 


(26) 


and the singlet/triplet parts of the pair potential expan¬ 
sion coefficient matrix 



(27) 


we can finally cast Eq. (20) into a convenient form: 


(at 

9k ) t '( 

<*. 

£ $ 

^) f - 

_ £ t J-(g£ 

9 k ' ) 

— Ek^kk 

>Tz , 








(28) 

with 








at 

= (% f 

ii 

u k 


■ 1 v? V? 

2t 

V k ■ 

..) T 

(29) 

9k 

=( 

• ii 
v k 

* 2t* 

V •• 

■ 1 4 T * 

* u T 

...) T - 

(30) 


The superscript 1,2,... means we have ordered the Bloch 
vectors and bands in some way. The precise way of order¬ 
ing is unimportant. Note that the set of {g k } solves the 
eigenvalue equation similar to Eq. (21) with the negative 
eigenvalues —E k while the set {g k } corresponds to the 
eigenvectors with positive eigenvalues E k . The elements 
of the set {g^} are the SC KS orbitals of SpinSCDFT 
in the normal KS orbital basis. We may easily repre¬ 
sent the densities using the normal state KS orbital basis 
{^icr(f)} for example 


n{r)= Y <Pi{ra)(n i:j ) C r (T np j (ra') , (31) 

icrjcr' 


and similar for m(r) and x{ r i r ') where x(r,r') is ex¬ 
panded in (p*{rcr ) and 0*(r'a'). The coefficients read 


(nij)aa' = (cr 0 ) aal Y(( ul k)* u k r fp( E k) + 

k 

+4>f')* M-E k )), 

(■ mij)av' = (S) CTff < Y(( U k)* U k fp( E k) + 

k 

+ v% k( v k T ) fp{~ E k)) , 

(x« w = (*w + 

k 

+<(<')* fp(-E k )) . 


(32) 


(33) 


(34) 
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We want to stress that we have not performed any ap¬ 
proximations so far and the SC KS system reproduces the 
exact interacting densities of the Hamiltonian of Eq. (1). 

b. Singlet Superconductivity Due to the antisym¬ 
metric structure of the fermionic wavefunction and the ef¬ 
fectively attractive interaction, in absence of magnetism, 
the singlet solution always leads to a more stable SC 
state. Known SC that feature a triplet pairing all share a 
very low critical temperature less than a few Kelvin 9-15 . 
In presence of magnetism, as we have seen, the spin is 
not a good quantum number and singlet/triplet compo¬ 
nents mix. Since the triplet pairing channel seems to be 
rather unimportant for many systems, it is of use to de¬ 
fine a singlet approximation, in which it is completely 
disregarded. 

We therefore make the assumption that our pairing 
potential has only the singlet component (marked as a 
subscript S in the KS potential). In addition, we assume 
a collinear spin structure in the normal state part of the 
Hamiltonian: 


A? A B 


$iZ\ s s 


ccrcr 

°ij 


°ij u o 


(35) 


Then, we observe that spin becomes a good quantum 
number in the SC KS system. This follows because the 
KS Hamiltonian matrix elements can be brought to a 
Block diagonal structure in Nambu and spin space with 
two kind of eigenfunctions to each individual block. Con¬ 
sequently we re-label the eigenvectors with k —> k, p 
where the size of the set of k is reduced to half. Each 
block p is diagonalized as 


9kn 9k,-ii ) 


t ( sign(^)A= 

Vsign {p)A s ^ -£~^ T 


= S kk , ( 


kfi 


o 


with 


9k>, 


9 k ',-L 


0 

^k/j. J 


1/4 

i 1 —/4 


u kfj, 

^ kfi ' ' ' Tfc/x 


•it* 

% ' • ' 

kf 


■■■? 


(36) 

(37) 

(38) 


E k is an eigenvalue that may or may not be positive. 
However, we have introduced the SC KS particles in 
Eq. (13) with a positive excitation energy Ek^ so this 
fact requires further commenting. In the present situa¬ 
tion where the matrix elements of the SC KS Hamiltonian 
are block diagonal in Nambu and spin space we can show 
that if g k has the eigenvalue E k the “negative” labeled 
eigenfunction g k _ /x has the eigenvalue —E k _ A7 . Thus 
we still have the original redundancy in the eigenvalue 
spectrum but not in the same spin channel p. Instead 


E 


± = —E t 

kfi k,—fi 


(39) 


We conclude that to every k we have 4 eigenvalues of 
which 2 are positive. These positive eigenvalues are iden¬ 
tified with Ekfj ,. In the next Subsection IIIA2c after in¬ 
troducing the Decoupling approximation we will be able 


to compute these eigenvalues explicitly, and continue this 
discussion. 

c. The Spin Decoupling Approximation It is desir¬ 
able to reduce the effort to solve the KSBdG Eq. (36) 
further. A substantial simplification is the Decou¬ 
pling approximation 33,35 (or Anderson approximation 48 ). 
There, one considers only singlet SC and pairing between 
a quasi particle state ( ia ) and its time reversed hole state 
(— i, —a). Furthermore it is assumed that the basis {(pi a } 
approximates the true non SC quasi particle structure 
well enough. In the language of the our KSBdG Eq. (28) 
this reads 


£ij ~ £ ia&oo'&ij > ($ • A ~ ^1 A s ■ 

( 40 ) 

This type of approximation is inherent in the Eliash- 
berg equations as well as SCDFT functionals. It is also 
straightforward to include a diagonal correction 7 . In 
the form presented here we will call it Spin Decoupling 
Approximation (SDA). For each k and p, Eq. (28) re¬ 
duces to the 2x2 equation 


l v k k :„ 


( < 

■ l v~ k \ 


( ^kcr 

\ v sign((r)Z\g 


V-ko 
,~k* I 

l -k-o 


* sign(<r)A“ fci _ fe 

-k,k — £ -k-a 


(f 



(41) 


Here we have introduced a single spin notation v,^ a = 
Vu_ a and u k ° = u ka . The spin label on the coefficients of 
the Bogoliubov transformation always refers to the nor¬ 
mal state KS basis spin label and thus we use the spin 
notation p —► a. Note however that the spin label cannot 
be strictly identified with the spin of a SC KS particle. 
We will come back to this point later. From now on 
we we will use the notation A s sk = A B sk _ k = A s s _ k k . 
We may compute the two eigenvalues and eigenvectors 
analytically. From the high energy limit £ka + £ -k,-a ^ 
£k a —£-k,-cr w e identify the name ± for the two branches. 
The eigenvalues are 

Ei. = ( 42 ) 

El = 524=1=^1 V C-+;— )W- (43) 

In the spin degenerate limit, the + branch has always pos¬ 
itive eigenvalues E ka and it is clear which of the eigenvec¬ 
tors belongs to the first column of the Bogoliubov Valatin 
transformation. In the spin polarized case the situation 
is more complicated. Again, because E k(J = —E^ k _ a 
two of the four Bogoliubov eigenvalues to a given k are 
positive but without knowledge of £k a and A B sk one can 
not tell in advance which ones these are. The general 
situation is sketched in Fig. 1 for a constant A B sk and 
homogeneously splitting free electron gas. In the next 
paragraph we give a more detailed discussion of the Bo¬ 
goliubov eigenvalues E ka . 
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k 

a) 



k 


b) 


Figure 1. (color online) Sketch of the Bogoliubov eigenvalues E^ a for a free electron gas with a homogenous splitting eta = 
k 2 + sign(<j)p B Bo- We choose a constant A\ k > g B B a in a) and Al k = 0 in b). We plot the + Bogoliubov branch in red and 
orange for'fand J, and the - branch in light blue and dark blue for band 4, respectively. We indicate the Ektr in a) as thin dashed 
lines. In a), the + branches are strictly larger than the Fermi Energy Ef and thus constitute the SC KS particle excitations. 
On the other hand for Al k < g B B a as in b), the + and — branch partly swap their order. When E k1 . > Ef the SC KS particle 
excitations are from the — branch also. 


d. Eigenvalues in the SDA Our first concern is how 
to interpret the spin quantum number a of E k(J in con¬ 
nection with the underlying normal states e k a- 

First, consider the non-SC limit where 

A B sk = 0 : 2E ka = £ k(J - £-k-a ± I Eka + £-k-c r| -(44) 

This situation is plotted in Fig. 1 b). Note that if £ ka + 
£-k-a > 0 , than E ka = -£_ k _ a and if e ka + £- k -a < 0 
we conclude E ka = £ kcr - 

Second, consider the following case that occurs at any 
k 0 where £ ko -\ + £- ko i = 0. Given that we have an energy 
splitting £fc 0 -f — £- ko i > 2|Z\g fc | we find that both E_ ko ^ 
are negative. This means that according to the defini¬ 
tion in Eq. (13) to take the positive eigenvalues, both 
KS particles are from the a =f branch. It is not possi¬ 
ble to construct the Bogoliubov transformations in this 
case and in any case the 7 ^ state cannot be occupied 
twice. It is, however, possible to give up the requirement 
that all SC KS particles are positive and simply always 
take the + branch. Then we can say that 7 ^ creates 
a negative energy excitation which will be occupied in 
the ground state. By analogy with BCS, 7 I creates an 
electron like single particle state on the SC vacuum, this 
leads to the interpretation that, in the ground state, this 
k space region is occupied by unpaired electrons. A sim¬ 
ilar discussion can be found (still in the context of BCS 
theory) in the work of Sarma 49 . Similar to Eq. (13) we 
can redefine electron to hole operators at the price of 
changing the ground state energy. Because the ground 
state energy, in turn, cancels from the thermal averages, 
the expectation values computed with this theory do not 
depend on this interpretation. We want to point out that 
this discussion only applies when the splitting is larger 
than the pair potential. 

e. Eigenvectors in the SDA Furthermore we can an¬ 
alytically compute the normalized eigenvectors g k to the 


eigenvalues E kfi (a = ±). We introduce the notation 


9ka 


ka 
U ka 
-ka 
h—cr 


(45) 


to label the components which are given in terms of the 
eigenvalues and components of the matrix by 


v 


— ka 
k—o 


u ka 

u ka 


/ I £k °'\ 

\j\Et + Et k -X 

sign(cr) A‘ sk I \£-k-g + E%J 
sign (a) \E% a + EP k _ a \ ' 


(46) 

(47) 


Starting from a converged zero temperature normal state 
calculation, within the SDA the only remaining variable 
is thus the matrix elements of the pair potential A a sk be¬ 
cause the SC KS wavefunctions as well as the Bogoliubov 
eigenvalues are explicitly given in terms of it. 

It is important to point out that within the SDA A a sk 
can be chosen to be real 40,41 . This can be proved by 
exploiting the gauge symmetry of Eq. (41) under rotation 
about the r- axis. If the rotation is applied with a k 
dependent angle 6 k of 


0 k 


= arctan 


( 


3^ 

KAlJ 


(48) 


we get: 


£ k a sign(cr)Z\” fc 
sign(cr)A=* -£-k-v 

£ka _ sign (a)Al k 
sign(cr)Z\| fe - £-k-a 


(49) 


where A B sk = sign(3f?Z\| fc )|| £ R. Thus the (k,—k) 
matrix elements of our general complex decoupled pair 
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potential are gauge equivalent to purely real ones. We 
still keep a general complex notation for A s sk first, to 
investigate explicitly if self-energy corrections affect this 
conclusion and, second, to make it easier to extent the 
formalism to the case where the gauge symmetry does 
not have enough freedom to make all matrix elements 
real. 


3. Competition between SC and Magnetism in the SDA 

The SDA, as introduced so far, assumes that we com¬ 
pute SC on top of a (magnetic) quasi particle structure. 
Thus, for example, it does not allow magnetism to be 
suppressed when a weakly magnetic system becomes SC. 
In conventional SCDFT 33 ’ 35 this type of feedbacks can 
be safely neglected because SC changes the dispersion 
only for states very close to the Fermi level. The effect 
on the electronic density is thus negligible and so is the 
change in the normal state xc potential. However, since 
the contributions to m(r) are in general more localized 
at the Fermi level, assuming quasi particle energies £i a 
to be unaffected when SC sets in may not be reasonable 
for magnetic systems. 

We want to point out in here that it is also possible to 
keep the simple form of the SDA and include competition 
of SC and magnetism at the same time, by means of the 
following iterative scheme: 

1. Take the normal KS states {c pi a } and eigenvalues 
£ka as starting orbitals. 

2. Solve the KS-BdG equations in the SDA 

3. Recompute the densities n(r) and m(r) according 
to the Eqs. (32) and (33) 

4. Re-diagonalize the normal state KS equations with 
the updated densities (in particular changes in 
m(r ) may be of relevance) 

5. iterate from point 2. until self consistence is 
reached 

This procedure changes the meaning of the SDA during 
the iteration because we are self consistently updating 
the orbitals {(pi a } it refers to. 

B. The Sham-Schliiter Equation of SpinSCDFT 

So far we have presented the structure of SpinSCDFT 
with the focus on the electronic SC KS system. However 
explicit functionals for the ire-pairing potential A s sk have 
not yet been discussed. The derivation of the approxi¬ 
mations for the xc-potentials generalizes one proposed by 
Marques 50 in SCDFT and uses the Sham-Schliiter equa¬ 
tion of SpinSCDFT. This equation is based on the obser¬ 
vation that the parts of the KS GF and the interaction 
GF that correspond to the densities must be equal. Using 


the Dyson equation for a SC in a magnetic field starting 
from the SC KS system as the formally non interaction 
one we can relate the xc-potentials to an approximation 
for the self energy. Here and in the next section we 
present a derivation of an xc-potential for SpinSCDFT 
that generalizes the ones of Marques 50 and Sanna and 
Gross 51 . 

We introduce the GF with the r ordering symbol T 
and the field operators in the Heisenberg picture 

G(rr, r'r') = —(T!?(rr) ® ^(r'r')). (50) 

The imaginary time ordering symbol in Nambu space T is 
defined to act on every of the (4x4) components individu¬ 
ally which can be achieved by transposing in Nambu-spin 
space 

T [r' t') = 0(t— r , )'F(rr)®>F^(r , r , ) 

—0(t'— T)(lf ^ ' f (r'r , )®>P ^ (rr)) Ts,, . (51) 

We define the equal time limit in the —1, —1 component 
different to the usual one (that we use in the 1,1 compo¬ 
nent). The equal time limit of the time ordering symbol 
should be defined to recover the density matrix operator 
but according to the usual rule where the creation oper¬ 
ator is taken infinitesimally before the annihilator would 
lead to the form in the —1,-1 component. From the 
equation of motion we derive the Dyson equation start¬ 
ing from the SC KS system as a formally non interacting 
system 

G(r,r',w n ) = G KS (r, r', w n ) + 

Jdr x Jdr[G KS (r,r h uj n ) ■ X' s (n,r / 1 ,a;„)-G(r , 1 ,r( w„), 

(52) 

with 

r s (r,r',w„) = £(r, r',w n )— v xc (r,r') . (53) 

Here £ is the irreducible Nambu self-energy, where the 
electronic Hartree diagram was subtracted, and v xc is the 
Nambu xc potential 

v xc (r,r ') = 

5(r — r')[crQV xc [r) — S-£? xc (r)) $ • Zl xc (r,r') \ 

— $ • A™*(r,r') —6(r —r')(CTov xc (r) —S*-i? xc (r)) J 

(54) 

The SC KS Greens function satisfies 

jc\ri(iuj n b{r - ri)T 0 cr 0 - ff K s(nn))-G KS (ri, r', w n ) 

= 5(r - r')r 0 cr 0 . (55) 

From the equation of motion we can compute the SC KS 
GF. Because by construction the SC KS GF yields the 
same densities as the interacting system we can cancel 
the respective parts of the GFs in the Dyson Eq. (52) 


that correspond to the densities. The result is the Sham- 
Schliiter equation 


J dri J dr i G KS (r,r i,w n )- 

■E B {r 1 ,r' 1 ,ui n ) ■ G(r[, r'u n ) 
G KS (r,ri,w n ) ■ 
■E B {r 1 ,r' 1 ,ui n ) ■ G(r' h r,uj n ) 


= 0 


J Oi, —a. 


= o. 


(56) 


(57) 


For convenience the self energy is decomposed in a 
phononic part E ph (u) n ) and a Coulomb part E c (oj n ) : 


S(uj n ) = £ ph ((jJ„) + S c (w n ). (58) 


S(u! n ) has a diagrammatic expansion in terms of 
G(w„) 24 and can be even viewed as part of a Hedin cycle 
for a SC including phononic and Coulomb interactions 52 . 
We do not consider vertex corrections, thus the Coulomb 
self energy part E c is the electronic GW diagram 



As an interesting extension we could include parts 
of the vertex corrections that lead to spin fluctuations. 
These, in the form of an effective spin interaction, are 
discussed by Essenberger et al , 53 and the extension to 
the present spin dependent formalism is straightforward. 
As compared to the polarization corrections of the same 
order, the phononic vertex corrections are negligible 54 . 

Moreover due to the quality of the phonon spectra one 
obtains with density functional perturbation theory 37,38 
we do not consider further diagrammatic electronic 
screening and treat the phononic interaction in the 
Hartree-Fock approximation 


^ph (^n) ~ 

It has been observed that computing the GW quasi 
particle band structure in a metal gives usually small 
corrections to the KS bands (compare Ref. 55 Fig. 2), 
also densities result to be almost identical. Thus, at 
least in the spin degenerate case, the GW corrections on 
a KS band structure of a metal are usually neglected. 
For convenience we use a similar assumption for the 
spin part. This way we can drop the Nambu diagonal 
v xc construction from the Sham-Schliiter equation. 
Representing G KS (r, r', u> n ) and G(r,r' ,w n ) in the 
same basis as the Bogoliubov-Valatin transformations, 
i.e. essentially the normal state KS orbitals ($'T^(r)} 
with the pure Nambu and spin spinor wavefunctions 



<F KS (r) 
^ 1006 \’ ) 


^a,l ^Pia (^) \ 

5 a,-l V* ia {r) ) ' 


(61) 


Sorting the expansion coefficients of G KS (r, r ', w„) = 

E«a , ™'.yCjaV , k)C( T ') ® in similar 

Nambu and spin form we obtain the 4IV x 4 A" matrix 
equation 

-^G KS (w n ) • ^ ^ S ) t Q j • G(w n ) 

_ \ ( ( 0 ^ 1 { UJ n) 

~ ( n) ‘ o 

+^ P h(w n )^ • G(w„) (62) 

that we need to solve for $ • /X s . From here on we use 
$ ■ Z\ a and 3> ■ zX xo synonymously, i.e. the external pair 
potential is assumed to be infinitesimal. 

In the next section we reduce the problem to the sin¬ 
glet case and employ the SDA. Because we can solve 
the KSBdG equations analytically we obtain a potential 
functional theory and arrive at a functional form that 
is formally similar to the BCS gap equation. We stress 
that the methods presented here and in the next sec¬ 
tion could also be applied without the restriction to the 
SDA. However in that case the equations would have an 
implicit form and require a numerical solution of the KS¬ 
BdG equations. Such a general form would be of impor¬ 
tance in considering triplet superconductivity or to ac¬ 
count for pairings beyond the usual one of time reversed 
states (as would be needed for example to describe the 
FFLO state 5,6 ). A further discussion can be found in 
Ref. 56. 


C. Derivation xc— Potential 

The Sham-Schliiter Eq. (62) involves the interacting 
GF which is usually only available after solving the Dyson 
equation. In an approximate scheme this step can be 
avoided. The straightforward way is to replace the ma¬ 
trix G(oj n ) with G KS (w n ) on all occurrences. As was real¬ 
ized before 33 this violates Migdal’s theorem because there 
the vertex is compared with the polarization diagram 
of the same order. Thus the phonon vertex corrections 
are only negligible as compared to the Hartree exchange 
diagram with the full GF. To circumvent this problem 
some of the Authors introduced a procedure to con¬ 
struct a self-energy that does satisfy Migdal’s theorem 51 . 
Starting from an electron gas model with a phononic 
Hartree exchange diagram, this leads to excellent agree¬ 
ment with experiment while still retaining the numer¬ 
ically simple form of the Sham-Schliiter equation that 
is independent on G(w n ) and involves only Matsubara 
sums that can be evaluated analytically. The self-energy 
E ks (uj„) with G(co n ) replaced by G^ s (w n ) is the basis of 
all further improvements. In this work, however, we will 
not investigate the parametrization procedure. We will 
limit the complexity of the derivation by using assum¬ 
ing E(uj n ) ~ E KS (cj n ), where in J7 KS (w n ) the G(w n ) is 
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replaced by Gf?(u> n ). This will give inaccurate critical 
temperatures but qualitatively correct results. Thus we 
are left to solve the equation: 


“ n x ' 

- I Vr KS fw i (( 0 E^KA 

[ nh ^V r ? s_I,I K) 0 ) 

+EhK)) -G KS (u;„). (63) 


In this form the matrix elements of the SC KS GF in the 
normal state KS basis are given by 


Gff^n) = 


“ i- E k 


k \ v{ ® uf 4 0 V 3 k * 


+E 


1 

iw n + Kfe 



+ 

(64) 


We use = ( A A) T with the expansion coefficients 
A of ■Ufe(r) in <p ia ( r) given in Eq. (23). Similar for FJ?,. 
Further we assume the SDA for the rest of this paper. 
Results beyond the SDA are discussed in the PhD thesis 
Ref. 56. In the SDA the SC KS GF simplifies to 


E/K) = 


( K?lA 


E 


\u n — E?. 


o 


0 


“it K 


"Si.-i 


K7I 2 <5ij 


1 UJn — E^ 


lu; 


Icd^— £ 7 ? 

K?lE 




1 uj n — E? 

0 

0 


lcc^+iiA 


1 u) n -\-E? 


(65) 


This form and any further formula based on it use the 
components of the SC KS wavefunction as given in the 
Eqs. (46) and (47). In the Dyson equation G -1 = 
G KS — E we see that we need to compare the self¬ 
energy contributions with the inverse SC KS GF. Invert¬ 
ing Gfj{ijj n ) with obtain 

(G KS )^ 1 (w„) = Sij (iu) n T 0 cro - f 4 + 2 '~'‘ i ) T z<r o 

~ C lt 2 £ ~‘ l )'r z a z ') +6 i ,- j ((iT y )(ia y )$iA B si 

+T x {ia v )iAA B si ) . ( 66 ) 

Here we see that self-energy contributions oc T z a o change 
the average spin Fermi level £lt+ 2 £ ~‘ 4 ' = 0. Similarly con¬ 
tributions oc r z a z change the splitting of single particle 
levels. It has to be understood that these are global prop¬ 
erties of the band structure, meaning that the full £; CT 
dispersion has to be integrated to obtain N electrons per 
unit cell. If the interaction changes dispersion and occu¬ 
pations far away from the Fermi level this may still cause 


a shift of the original Fermi level. An clear cut exam¬ 
ple is the following: In the context of SC one often em¬ 
ploys the Eliashberg function o?F(f2) which is the Fermi- 
surface average of the electron-phonon interaction 25,27,28 , 
to describe the electron phonon interaction. This func¬ 
tion is assumed to apply equally to all states, also those 
away from the Fermi level. This is a good approxima¬ 
tion only if corrections of the Fermi level are excluded a 
priori (electron-hole symmetry), otherwise under this as¬ 
sumption the correction to the Fermi level Eit 2 e ~ 11 and 
the splitting ait ~ 2 g ~ u would show a logaritmic divergece. 
As commonly done in Eliashberg theory, where the same 
effect occurs, one then excludes self-energy contributions 
oc t z . We will assume the same approximation. As the 
Hartree diagram is proportional to t z is thus not consid¬ 
ered. While the expected Fermi energy shift is negligible, 
corrections to the spin splitting £lf ~ 2 S ~ U cou i ( ] kg 0 f rele¬ 
vance. However due to the extreme additional numerical 
complexity of considering the true full electronic state de¬ 
pendence of the electron phonon interaction we leave this 
to a future project. We compute the self-energy matrix 
elements in the SDA from the Eq. (60) 


qka 


X \ u ka\ 2 M ph (f2 q , E krT , UJ n ) 

(67) 

EhicrjcrE") = E, 9jkcr9i.-k.-cr X 

qka. 


X Wfe“ {Vk-v)*Mp» (fig > E k *, w„) 

(68) 

Eh x 

qka 


X ( u ka)* V k-aMph(G q , —Ekcr, W n ) 

(69) 

Eh iaja' A) = ^ aa ' ’ S ^ J 9kia9jkc X 
aka 


x\4:\ 2 M ph (n q -E^,co n ). 

(70) 


From the hermiticity of H^ s ph of Eq. (10) comes g°L q {r) = 
(g“(r))* and thus the electron phonon interaction matrix 
elements 


9ija = [ dr E E( r ) ‘ ‘ EE)<(r) ( 71 ) 

J - n 


a=0,z 


have the property g\^ a = gA a . Moreover g\- a oc 
£>ki,kj+q which is expected from the lattice translational 
symmetry 37 . The Matsubara summation M ph (I2, E,ui n ) 
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is evaluated with the result 

Y ph (f2, E,ui n ) = -qY,m.-e\ 


(3 *-f iu/ n > - E i(w„ - w n >) + J7 

n' 

np{f2) + f 0 (E) 


(72) 


M ph (fl , E , ui n ) — 


£2 — E H - 

npjfi) + fp(E) + fp(E)+np(-f2) 
£2 — E + iCJ 77 , £2 E — i ujn 


(73) 


(74) 

= y ph (i?,£;,a; n )-y^(i?,-£,w n ) ) ( 75 ) 

where fp(E) and np{D) are Fermi and Bose functions, re¬ 
spectively. The Coulomb self energy parts on the Nambu 
off diagonal with the diagram of Eq. (59) are 

^ S ViA“n) = — 6<r,-(r' E W lkl-k a ^ >< 

ka 

*4>k-TU(E^), (76) 

vr&Aun) = -**,-*> E x 

ka. 

x4TVk-:M-ED- (77) 

with the static screened Coulomb matrix elements 


TT^stat 

kik^k^k^fjfj 


' = J dr ' ^ CT ( r ) x 

e - 1 (r, r', 0 ) 


r — r 


^k 3 Ar')-^'(r'), (78) 


with the inverse dielectric function e _ 1 (r, r', 0 ) that 
is accessible in many electronic structure codes 57,58 . 
e” 1 (r,r',0) is often calculated within the RPA which 
yields very good results for metals in general. As we 
have pointed out, terms proportional to r 2 i.e. contribu¬ 
tions (E«f — E^ ’ x ) are dropped from the functional 
construction. 

Because of the gauge symmetry discussed in 
Sec. Ill A 2 c, we expect the equations for A B sk and A B sk to 
be similar. Thus we proceed and evaluate only the 1, —1 
component of the Sham-Schliiter equation (62) in SDA 
and arrive at 


M k k] Z k k Al k + M% k k AH = © fc ,_ fc + . (79) 

Here T>k,-k are the purely phononic contributions due 
to the Nambu diagonal self energy parts ro(A' [ ') t f 1 ’ 1 + 
E K ^ 1 L )- £fc,-fc is due to the Nambu-off-diagonal self 
energy contributions and contains the phononic interac¬ 
tion along with the Coulomb potential on the same foot¬ 
ing. The coefficients 


M, 


k-k 
k, — k 


M 


fk,—k 

-k,k 


3 E < 8 °) 

^ na 

(E s KA->»> <3 !(81) 

^ na 


are the Matsubara summed SC KS GF parts. Note that 
M ,k ! k k A e sk oc \A e sk \ 2 Al k so the Sham-Schliiter equation 


in the SDA is unaffected by the phase of A* sk , as expected 
from the gauge symmetry. 

T>kk’ and C^k 1 have non vanishing matrix elements 
apart from k, —k. These are not included in the SDA. 
Other SC theories such as Eliashberg and spin degener¬ 
ate SCDFT are build on similar approximations and from 
the quality of the results one obtains, we conclude that 
such corrections are in general not important. 

Another interesting aspect of the functional construc¬ 
tion to observe is that a self-energy part showing tx 
triplet symmetry appears, that means the spin inverted 
Nambu off diagonal components are not equal and of op¬ 
posite sign 


E-^ + E^^O. (82) 

These self-energy part leads to non-vanishing functional 
contributions in £k,-k in the singlet channel. We call 
these contribution intermediate triplet contributions. We 
have investigated the effect of removing them and found 
that this has essentially no consequence in the numerical 
calculation for a spin independent coupling (see part II 
). In addition we note that similar to the matrix elements 
k' ^ —k, the diagrams generate triplet contributions that 
cannot be incorporated into the SDA. This also means 
that the terms 

^(^G KS1,1 K)-F KS1, - 1 ( Wn ). 

a n 

■G’ KS “ 1 ’“ 1 K)b,_ <r ^ 0 (83) 

E(a E^ 31 ’" 1 ^)-^ 3 " 1 ' 1 ^) 

a n 

•G KS1 ’- 1 K)) CTi _ ff ^0 (84) 

are not zero as, on the other hand, one would expect for 
a singlet SC. This fact simply means that ignoring the 
triplet components from the external potential is not con¬ 
sistent, in presence of a magnetic field, because a triplet- 
singlet coupling exists at the level of the roc-potential. As 
discussed earlier (Sec. Ill A 2 b) , it is not clear in which 
cases triplet effects become relevant. However, since ex¬ 
perimentally triplet SC is only observed at very low tem¬ 
perature, in high temperature regimes disregarding all 
triplet components should be safe, we will show in II 
when we investigate the influence of intermediate triplet 
contributions, at least, that they are small. 

Within the SDA the SC KS wavefunction components 
v k-a > u ka are explicit functionals of the potential A B S . 
Thus, left and right hand side of the Sham-Schliiter equa¬ 
tion (79) are equally non-linear functionals of the poten¬ 
tial A®. We interpret the Sham-Schliiter condition (79) 
as 


S p [A B s ]-Al = 0 Sf) = S% + S* hfi + S* 0 + S% . (85) 

Here Sf-A B s is equivalent to -~{A B sk M k ’~ k + M ,k j k k k A B s * k ), 
Sp -A b s = D k ,-k and ■ A e s = <t k ,-k • The non- 
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linear Sham-Schliiter operator contributions are given by 


CM £ \ ' ( i^ka ~r t' — k—a) D t ^ 

S P kk' \ F +_ F ~ 12 ^ E ka^ka) 

a ^ial 


£ — k — a 


+2\u[ 


fc+l 2 U ,- fe+l 2^p s(jB a ff)jB o 


ka I 


u k-o 


( 86 ) 


and 

S 


1 


P kk' ^ > \E+ — E~ 

qk2cr0t.\Qt.‘20i-3 ' ka ko 


sign(a 3 ) 


: -k 2 OL2\2\ -ka. 1 \2\ q |2 

^k 2 — cr I wfc —a I \9—k2, — k, — cr\ 


x ((|t 

+ (l«fa-TI 2 |«fc-“ 1 | 2 |5l fc> - fca> -J 2 + 


(87) 


The term S? due to the Nambu diagonal acts to re¬ 
duce the critical temperature. In the Refs. 33 and 35 this 
term was scaled down by a factor of 1 /2 in the functional 
construction to compensate for a systematic underesti¬ 
mation as compared to the Eliashberg critical tempera¬ 
ture in the phonon only case. In Ref. 51 a SCDFT func¬ 
tional is constructed, by using a proper interacting GF in 
the exchange self-energy of Eq. (60) , therefore removing 
the necessity to reduce the repulsive Sp . Having in 
mind to generalize this functional to SpinSCDFT, in the 
present work we decided not to use the scale factor. In 
part II we find further indications that this scaling may 
also effect the robustness of the SC state against a mag¬ 
netic splitting. The predicted critical temperature will be 
too low as compared to experiment but the correctness of 
the qualitative behavior of the theory will be preserved. 
The Nambu off-diagonal contributions that derives from 
the phonon interaction then reads 


CE 

p h Pkk 


' = -£ £ 




qa o:iQ2Q:3 


q —q 

9kk'a9 —k, — k' , — a 

sign (a 2 ) ^ — E^i J 


kai | 2 |^—fca 3 |2 


U ka U k-c r U ka U k-a A 


A; 


s k‘ 

U 

A sk‘ 


-^f2 q ,E^,E^,E«*).. 


( 88 ) 


and the contribution that derives from the static 
Coulomb interaction reads 

y^ _ sign(a 2 ) 

a 0,0203 


CE 

A C0kk' 


\E±„ - E7. 


x(w&,-k,-k> a _j4?\ 2 K-T\ 2 


+u 


kaii 0 ,—kai* ni ka3 o} —ka3* ^sk' ^jy-stat 


ka ^k—a 


U 


ka V k—a 


A 


kk'—k. — k' 


sk' 


Xic^W)- 


(89) 


The functions P s , L and L c coming from analytic Mat- 
subara summations, are given in the Appendix A, to¬ 
gether with a discussion on some limiting cases. 


1. Description of the Second Order Phase Transition 

If the SC transition to the normal state is of second 
order, x( r , T ’ / ) is assumed to go to zero continuously 
upon approaching the critical temperature. From earlier 
work 49 in the BCS framework, we expect this to be the 
case in the low magnetic field part of the phase diagram. 
The formalism in the SDA is built on the potential A® 
not the order parameter %. We thus need to proof that 
a second order phase transition implies also a continuous 
vanishing of the potential A®. We note that in the SDA 
it is sufficient to show that the expansion coefficients of 
X and A s s in our normal state basis are of the form 

Xska,-k,-a = a k~ aA lk - ( 90 ) 

where ct k ' a is some function of A® and show that 
a k’~ a ( A s) 7^ 0- Given that this is the case, in 
the limit |A®| —> 0 only linear order terms in the Sham- 
Schliiter equation are relevant. Then, at a second order 
phase transition T c can be computed from the condition 
that the matrix limi^si^g 08 [A®] is singular. 

Coming back to Eq. (90) and using the SDA together 
with Eq. (34) we see 


fp(Ka) ~ fp( E ka) 
I Et a -E^\ 


(91) 


Clearly, at T > 0 a k can only be zero if E kcr —E k(J —► 0. 
Taking the respective limit 


1 - cr, —cr 

hm a k ’ 

E+ -E7 ->0 

ka k a 


_ ( 92 ) 

2 cosh(/3(£J^ 0 . + E ka )/2) 

< 0 (93) 


which is the desired result. We may thus use |A®| —> 
0 instead of y s —► 0 at the point of a second order 
phase transition. We sketch the function a k ’~ a using 
A = 0(25+ + 25*")/2 and B = 0(25+ - E ka )/2 in Fig. 2 
Note that while aj: a is strictly non-zero if @{E k + 
E k )/2 1 then a k '~ a is exponentially small in the 

range \B\ <C |A|. . We thus observe that the order pa¬ 
rameter Xs{ r i r ') i s only weakly dependent on the poten¬ 
tial matrix elements A®, that correspond to states below 
the splitting energy A. Still, this does not invalidate the 
conclusion that at any finite temperature a continuously 
vanishing order parameter implies a continuously vanish¬ 
ing pair potential. We thus expect that (at low splitting) 
we can use the linearized Sham-Schliiter equation (85). 
In the following, we use a breve on top of linearized enti¬ 
ties such as Sp c = lim^i^o 0g[A®] and Eq. (85) can be 
solved from the condition 


det Sp c = 0 , (94) 

where 0 C = 1/T C . The right eigenvector of Sp c is propor¬ 
tional to A®. To compute the small A® limit of Sp c we 
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Figure 2. Sketch of the function /3a^’ _<T (Zia) multiplying 
Al k into the coefficients for Xs{ r , r ')- While being always 
non-zero, the coefficient function at B = 0 behaves as 
— (/3/2)/cosh(/3(i5^" 0 . + E kcT )/ 2) and thus becomes exponen¬ 
tially small with decreasing temperature. 


first investigate the behavior of |u£“| , a: 

separately where we find 

l/iaj^o = ^o:,sign(Efc,,+£_*_„.) ) 

| Jj™ Q I V k-a I = ^c«,-sign(E*.„ +£_*_„) > 

|J™ 0 ' E7fcOr _ 

— &a,-s ign(e kc7 +e- k -A £ - k ’- a ' 

Also we see that 

lim \u k £\ 2 \v7 k *\ 2 = lim u k k %vT k f = 0. 

s straightforward to arrive at 

^0 kk 1 = ~ ^^'PsiE-k^i — £ —M-) 


Thus it is 




and 

^ fefc' “ 

- E (ifffclt I 2 Efct. £ lt. £ fct) 

fe t i " fc d 

+L(f2 q , £ k f, — £q, £fc-f) — L(n q ,£ k f, £q, 

—L(f2 q , Efc-f., — £q, — E-fcq.)) + 

+ |s , !.fc i _4| 2 (-£'(^qj e_ fe 4., e_4, £-*4) + 

-(- Z/(12g, £—/4, £—4? £—A4) £—^4, £ 

-L(J? g ,e_fe4.,-e_4,-efcf)) . 

Moreover 

Off 

pi>0kk' 

-*e i £ ^t, 4v - 

^|£fc't + £-fc'4l 

+L(f2 q , £-ki, —£fcf)) 

and 

_ _ 9 TT^stat 

^ ** kk', — k,—k '^ * 

j —e - 


, | ,4, 


,2. Non-linear Gap-Equation 

Far from T c or in those parts of the phase diagram 
where the SC transition is of first order we need to use 
the non-linear Sham-Schliiter equation, because a solu¬ 
tion with small |Z\® fe | may not exist. The most common 
method to solve an equation of type Eq. (85) is to use 
an invertible splitting matrix S and cast Sp [A®] • A B S = 0 
into a fixed point problem 

A® =JC s [Al\-Al JCs = S~ x ■ (Sp + S) . (103) 

This is the gap equation of SpinSCDFT. In the spin de¬ 
generate limit the choice S = —S “ leads to the SCDFT 
gap equation given in Ref. 33. We point out that while 
we can show that all So, , <0 at £ ka + £-k -a = 0 

P k,k 

S?k k ~ ex P {—\P{ £ ko — £—fc,—o-)) and is thus a numer¬ 
ically problematic object. In the implementation that 
we describe in detail in II we find that a good choice is 
S = —S™(£ka = £-k,-a )• Obviously in the spin degen¬ 
erate limit we recover the formulas given in Ref. 33. In 
part II we will also discuss the properties of the splitting 
versus temperature diagram for a simple system in detail. 


IV. ELIASHBERG EQUATIONS 

In the KS-SpinSCDFT formalism, interaction effects 
are mimicked by the ccc-potential that is an (implicit) 
functional of the densities. While the functional con¬ 
struction and the additional complications of the SC KS 
system pose additional algebraic complexity, the result 
is a numerically cheaper computational scheme. This is 
owed to the fact that Matzubara summations in the self¬ 
energy are not computed numerically but absorbed into 
the analytic structure of the a;c-potential. Likely, the 
knowledge of the interacting self-energy is essential to a 
future improvement of the presented functional. The self¬ 
energy Eq. (53) in turn is constructed via diagrammatic 
perturbation theory using the electronic and phononic 
GF similar to Sec. IIIC, and involving the solution of 
a Dyson equation. In the present section, we develop 
this direct many-body scheme to obtain the electronic 
GF. The final set of equations generalize the ones of 
Eliashberg 28 and we refer to them with the same name. 
Ref. 24 discusses similar equations in a different notation 
with a limitation to isotropic system with a homogenous 
splitting parameter. 


A. Solving the Dyson Equation 

The starting point of the derivation of the Eliashberg 
equations is the Dyson equation of a SC Eq. (52). We 
represent it in the basis of normal state, zero temperature 
KS orbitals {<Pia(r)} defined in Eq. (22). We use the 
Nambu-Anderson 40,41 notation similar to that used in the 
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functional derivation and in Eq. 62. The Dyson equation 
reads 

K) = Gf/K) +Ypf k s (co n ) • • Gij(uj n ), 

kl 

(104) 

where Gff is the SC ICS GF and 17^(w„) = £ iclj (u„) - 
v xc ij where S xo ij(co n ) is the Nambu exchange and corre¬ 
lation self-energy that also includes the phononic Hartree 
diagram 52 . u xc jj are the matrix elements of the xc- 
potential of the SC ICS system. Note that the SC KS 
GF is not diagonal in the space of Similar to 

our approach in SpinSCDFT of Section III we assume 
that {0 ia (r)} is a good approximation to the quasi par¬ 
ticle state 59 , i.e. Eh(uj n ) and G i3 -(w n ) are essentially di¬ 
agonal. We use similar diagrams (Eq. (59) and (60)) as 
for the functional construction of SpinSCDFT in Subsec¬ 
tion III C namely the phononic and Coulomb exchange 
diagram. Again similarly (compare Sub. Ill C) we drop 
the Coulomb corrections on the Nambu diagonal that 
add to the xc potential. Further we assume, as in the 
SDA of Sec. Ill A 2 c, that the pairing occurs only be¬ 
tween time reversed states 48 . This means we only con¬ 
sider singlet SC. Starting from Eq. (104) in the form 
Gij{u; n ) = (Gf/ _1 (w„) - r"j(w„)) , under the men¬ 

tioned approximations, the Dyson equation is a 4 x 4 
matrix equation that can be solved analytically. Note 
that here we do not substitute the SC KS GF for the in¬ 
teraction GF in the self-energy (as was done in the func¬ 
tional construction of SpinSCDFT of Sec. IIIC). This is 
the main difference in the two approaches so far. 


1. Analytic Inversion of the Dyson Equation 


SE part 

G KS 1 part 

Basis vector 

Eliashberg 


icj n 

TOITO 

Z k (gJtl ) 

Afc (Wn) 

- 

T~0 O' z 

Afc(t^n) 

( W ») 

(^fct H - £—ki)/2‘ 

T z CFq 

£ k (yGJn ) 

E k ( w ra) 

(£k- f ~ E-fe4.)/2 

T z Vz 

J k (^n.) 

^ A k) 


(i Ty)(i<Jy) 

^*(w„) 

< 

i %A s sk 

T x (ia y ) 

4fK) 


Table I. Self-energy contributions, the variable of the inverse 
SC KS GF which they add to and the basis vector. E.g. along 
the to oo direction in Spin and Nambu space i uj n + E k (io n ). 
In the last column we give the related Eliashberg property. 
Note that A% ~ rJf A + iX)? A and Af*(w„) ~ Ef A - iXfA 


responding scalar self-energy components read 


( w n) — ^ a - 0 (^ra' w n) x 

aak'q n' 

X \9kk'<T^G } :, a k , a [y>n') 

A k^n) = ^ XJ 0 Y, X 


(106) 


tTOtk'q n' 

x l5fcfcvl 2 Gfc; CTi fc/ CT (w n ') 


(107) 


^ n y^(' r z)aa-P 0 .-0(^71' W n ) X 

aak'q n' 

X \9kk'cr\ 2 Gk'' < j,k l a( UJ n') 

S k(<^n) = ^ XJ £ Y X 


(108) 


aak'q n' 

x| 4J 2 ^C,^K')- (109) 


The easiest way to invert the right hand side of the 
Dyson equation 

Gij^n) = (Gff 1 K) - ^K))- 1 , (105) 


Note that A k (cu n ) stands out in the sense that the SC KS 
GF has no contribution along this direction in Nambu 
and spin space. On the Nambu-off-diagonal we similarly 
introduce 


is to identify contributions of the self-energy that add 
to a given variable of the inverse SC KS GF Gff 1 (w ra ) 
of Eq. (66). We summarize these self-energy contribu¬ 
tions in Table I . This means we decompose the Nambu 
and spin matrix E kl (u> n ) along the vectors r 0 <To, t z <7o 
and so on. Then, we name the self-energy contributions 
according to the property of the SC KS GF they add 
to in Eq. (105) and indicate the property in the super¬ 
script. For example the Matsubara frequency variable 
of the inverse SC KS GF points along the Toco axis in 
spin and Nambu space. Correspondingly the self-energy 
part along basis vector is referred to as E k (u> n ). In the 
following we use \g q kk ,J 2 = \g k , q kr7 \ 2 , D° _ q = D°_ qq and 
Wkk'k'kaa = Wk'kkk'aa > T h en the equations for the cor- 


?■(».) =-£ l£ * 

ak' n' q 

X 9kk'a9-k-k' -a + ^kk',-k,-k ' (7 _ cr ) X 
X y~l(' r a:)aaG fc ; a . _ k , _ a {w n ') (HO) 

a 

irw =-£ \ -«,) x 

ak' n' q 

9k'ka9-k’-k,-a + W k'k,-k’ -k a _ a ) X 
X ^~^(i r i i)aaG k ’ lcr _ k , _ a (u] n ') ■ (HI) 
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Here we introduce B k (uj n ) = S k A (uj n ) + iS k A (uj n ) and 

i^ A K) 

^K) = - E ^ E (E - w «) >< 

ah' n' q 

X 9kk'a9-k-k',-a + W kk'-k,-k' a _ (J ) X 

(112) 

SJK) = £ ± E ^(E - "J x 

crfc' X n' q 

X 9k'ka9-k' -k-a +W%t,-k'-k <7 _ ( ^ X 
X ^'k'a,-k',_a( Un ') ■ (H3) 


If both E® A and T, k A are real B k is the complex con¬ 
jugate of Bk . Further, for the same reasons discussed 
in Sec. Ill A 2 b, we do not consider the possibility that 
triplet self-energy contributions appear. It is important 
to remark that, just as in the usual derivation of the spin 
degenerate Eliashberg equations, the k dependence of all 
self-energy parts is generated via the k dependence of the 
Couplings \ g kk , a \ 2 and in addition W k , k _ k , _ k on the 
Nambu off diagonal. 

Introducing the mass renormalization function Z k (u} n ) 
as 


Zki^n) — 1 +i ljfc(Wn)/Wn, (H4) 

we can rewrite some of the above equations by including 
Z k {uj n ) into the self-energy parts: 


Z\f(w n ) = B k (uj n )/Z k (Lu n ) (115) 

Z\fK) = B* k {co n )/Z k (co n ) (116) 

£k(w n ) = ((£fef + £- k i)/2 + Z k (uj n ))/Z k (uj n ) (117) 

4K) = ((fffct - £ -k \)/2 + Hi(v n ))/Z k (u) n ) (118) 

A%{u n ) = A%(u n )/Z k (u n ). (119) 


Then by introducing the abbreviation 

5to(w„) = ((i k (u n ) +sign(cr)Afc(w n )) 2 

+A k ( w «)2lfc*(w ra )^ , (120) 


and suppressing the arguments ui n , we arrive at the for¬ 


mulas for non-vanishing SC GF components 

Q 1A _ 1 + a{e k + sign(g)A%) 

ka ’ ka 2^ ka Z k ^ iw n - sign(cr)Jfc - ot$ ka 

( 121 ) 

q-1,-1 1 y^ gfc- CT +a(4-sign(o-)A^) 

kaM 2$ k -<rZk^ / iuj n + sign(a)J k + a$ k - a . 

( 122 ) 

qi,-i 1 v- sign(g)a^g 

ka ~k~ a 2$ ka Z k 2-'i U ) n - sign (a)J k - a^ k(J 

(123) 

q- i.x 1 y- sign(o-)aZ\f 

ka-k-a 2$ k _ (T Z k ^i U } n + sign(a) J k + a$ k - a 

(124) 


We have thus expressed the GF in terms of the self-energy 
components (Eq. (115) to (119)) explicitly. The coupled 
set of equations Eq. (115) to (119) are the Eliashberg 
equations and have to be solved according to the scheme: 

1. Start with the coupling matrix elements g k i ka and 
Wk'k-k' -k and an initial guess for the self¬ 
energy components A k , A k *,i k , J k and A k . 

2. Evaluate Eq. (115) to (119). They are closed in 
the sense that inserting the equations of this sec¬ 
tion Af, A k *, i k . J k and A k only dependent on each 
other and the coupling matrix elements g k i ka and 

TJ/Bta* 

yv k'k _ k' _ k a Q- ’ 

3. Construct a new self-energy and iterate from point 
2, up to self-consistency. 

A k is a peculiar object because it generates a spin im¬ 
balance in the particle as compared to the hole chan¬ 
nel. To understand the effect of A k consider the follow¬ 
ing self-consistent cycle. We start the iteration of these 
equations with A k = 0 and E k A = 0. Then follows 
Gfc'U.-a = G- k - a _ ka which results in B* = B k and 
no self-energy part E k A is generated. Further, because 
G k aka = -a we then that ^k i s proportional 

to the difference of the interaction in the spin channels 
A k oc \g kk i^\ 2 — \g kk 'i\ 2 - If now the interaction is indepen¬ 
dent on the spin channel \g kk ,^\ 2 — \g kk '^\ 2 = 0 then A k 
also remains zero and we are at our starting point. Thus 
we conclude that for spin independent couplings S k & and 
A k remain zero during iteration. If the interaction is spin 
dependent \g kk ^\ 2 — \9kk'i\ 2 ^ 0 the self-consistency it¬ 
eration will generate a spin imbalance in the GF. This is 
not surprising because the up and down single particle 
spectrum is altered in a different way by the interaction. 
Then a non-vanishing E k A cannot be excluded. 

For future reference we extract the renormalized en¬ 
ergy dependence £ k of the GF as it appears in the self¬ 
energy Eq. (106) to (109) and Eq. (112) and (113). With 
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the abbreviation 

afc(w n ) = (A k ~) 2 + AlAf + uj 2 — ( Jk ) 2 (125) 

we obtain (6 = 0, z) 


^2( T b)aaG kaka — 

a. 

^(r b ) Q « a^4-sign(CT)(J fc +ai£)) 


=E- 


£ 


fc a fe^ c ;(iw„ +£^ 


2a 

sign(a)' 


(126) 


and 


_k_ a (Un) — 


= 4 e 




k a a k — a(iisj n Jk + A k i k ) + £ 


(127) 


B. Analytic Integration of the Energy 

In a numerical solution, the equations (115) to (119) 
have to be iterated until self-consistency is reached. Each 
self-consistent step requires to perform Matsubara sum¬ 
mations in addition to the k space summations which will 
be numerically demanding. 

Note however that the k space summations can be 
avoided using an approximation that is very common in 
the context of Eliashberg theory which is essentially to 
replace the couplings with their value at £fc(w ra ) = 0. 
The reason why this is sensible can be understood from 
the GF. From the above equation (126) one can easily 
see that G^ kj (u„) - G k ^ k ^(uj n ) behaves as (4k)) 
for large £fe(w n ). In turn, G^ fc£7 (a; n ) + G k ^(u} n ) and 
the Nambu off-diagonal parts G k ’~Z k _,k) behave as 

(4k)) 2 for large £fc(w„). Using this insight we 
see from the Eqs. (106), (107), (112) and (113) that 
4k), A k (uj n ), A k (wn) and A k *(uj. n ) are almost inde¬ 
pendent on the space k belonging to large e k because 
its contributions are suppressed by a factor (£fc(w„)) 
Thus these quantities can be computed replacing the cou¬ 
plings with their value at £fc(w n ) = 0. 

With the integrand behaving as (£fe(w n )) , the con¬ 

vergence of the Brillouin zone integrals in E k (uj n ) and 
E k (co n ) depend on the fc-dependence of the couplings in 
an essential way, even on k that correspond to a large 
£fc. In particular, in absence of any k dependence of 
the couplings E k (u n ) and E k (oj n ) diverge logarithmi¬ 
cally. From the physical point of view E k (uj n ) shifts 
the position of the Fermi energy and E k (uj n ) the mag¬ 
netic splitting of quasiparticle states due to many-body 
interactions. These terms are zero if the system shows 
particle-hole symmetry and small in general (see also the 
discussion in Sec. IIIC). Therefore we will discard these 
contributions completely and replace the couplings with 


their value at 4k) = 0 entirely, reducing the compu¬ 
tational costs significantly. 

Another very effective simplification of the formalism 
comes from assuming the system to be isotropic in k. 
This means that the couplings will depend on k only via 
the quasi particle energy £ ka . Here, we introduce the 
averaging operation on a generic function F k(J on equal 
center of energy and equal splitting surfaces according to 


F a (e, J ) = ha(e, J)F kr7 


1 


g(e,J) 


E 6 (- 


£sign(<r)fc t F ^ —sign(CT)fc4- 


x6(- 


- £-„ 


— J)F ka , 


(128) 
e) x 

(129) 


where the number of states on center of energy and 
splitting surfaces is given by q(e,J) = /fc CT (£, J) 1- The 
subscript indices “fccr” on I k(J (E,J) indicate the vari¬ 
ables that are averaged. Note that we invert the sign 
of k for the a =\, part which makes Iko{£, J)F ka = 
J)F- k - a . Now we define the analog of the 
Eliashberg function c?F(lo) 25,2 ' . We are going to keep 
the state dependence k for a little longer, and eventually 
take only those k such that £fc(w n ) = 0. On the Nambu 
diagonal it appears the coupling function 

o?F°(e,J,e',J',Q) = 

e(E\J')i k>rJ {E\j')i ka (E, (130) 

Q 


and on the Nambu off diagonal 
a 2 F{ £ , J, £', J\ Q) = q(e', J')4v(e / , J') x 

x -4ct(£) J) E] 9kk'a9-k,-k' ~ ^?) (131) 

Q 

C stat (£, J,e', J') = 

q(e\ JO/fevk, J')iUe, J)W$%l-k‘,-k a ,_ a • (132) 

Note that in the above equations (131) and (132), the left 
hand side does not depend on a because the averaging 
leads to the same result for a =f or a =),. The summa¬ 
tion over k! and q in the self-energy Eqs. (106) to (113) 
are then transformed to integrals over e', J' and fi re¬ 
spectively. However, if the couplings loose their center of 
energy dependence £, the following functions only depend 
on the Matsubara frequency ui n (that we now indicate as 
the index n) and the splitting: Z n (J), A“(J), A®(J) and 
A%*(J). With £fc(w n ) = e/Z u and J k (ui n ) = J/Z n we can 
compute analytically the integral over the center of en¬ 
ergy £ of Eq. (126). Because the integrand decays faster 
than e^ 1 for large £, we may compute the integral 


®WJ) = 


rE- 


a(e- sign(cr)(J + A%Z n )J -i u n Z„ 


i{J)Z n sign(cr) (l U n Z n J + A^Z n £) + £ 2 


( 133 ) 
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as the sum of residues in the upper complex half plane. 
Since it is not clear which of the four poles will be in 
the upper half we compute all residues. Adding those, 
we obtain the energy integral in Eq. (106) and Eq. (107) 
with 


&n,cr{J) 


as 


\J ~{Z n 2 A^A™ - (i uj n Z n - sign(fj) J) 2 ) 


(134) 


m na (j) = 

'iu> n Z n - sign(cr) J 

©n,cr 

,(\u n Z n - sign(u) J 

A &r 


+7T1 


n,<r 

. /i io n Z n + sign(cr) J 


(3 

w n ,—<t 


-1 

+ 1 

+ 1)0 


4^ Z 
sign (a) 
A“Z„ 


H a (-£r^- e - 






7, 

^ > 
^n.—a . 


sign (a) 


,fiu} n Z n + sign(cr) J 

©ra. —ct 


-10 


'(*<5 


A u 7 

'Z-’n. — r 


sign(cr) 


(135) 


Further, for Eqs. (112) and (113), we integrate Eq. (127) 
in center of energy e. We define 


Wn(J) =Ei 


de 


a n Z% - a(i Lu n Z n J + eA%) + e 2 


(136) 


that is evaluated to 

m n {J) = 7ri(6 ra , T - 1 0(3(-<Z„ - e n ,t)) 

-6„ it _1 0(S(-i“Z„ + © n , t )) 
+6„ i4 - 1 0(3(i"Z n -©„.;)) 
-S n4 - 1 0(9KZ n + ©„,+))) . (137) 


We obtain the Eliashberg equations similar to their usual, 
spin-degenerate form 26,60 , that only refer to the GF im¬ 
plicitly 


z n (J) = i + J')on nl<7 (j') 


(138) 


K{J) = 


1 




4 Z n (J]l P f+ sign(tr) 




A ^ J) ~ 2Z n (jJ dJ 'P 

n' 

xZ n ,(J')A%,{J’)% l '(J') 

K*{J) = -b7Tn/ dJ, AE £ «AW^) 


2 Z n (J)J P Z-f 

n' 

xZ n ,(J')A%(J')Vl n ,(J') 


(139) 


(140) 


(141) 


where 


a: 


U'W J, ) = fa 

,,n'{J,J') = JdQ 


217 o?F°(0, J, 0, J', 12) 
(w n - uv ) 2 + J7 2 
2l7a?E(0, J,0, J',17) 


(w„ - w„/) 2 + 17 2 
+C atat (0, J, 0, J'). 


(142) 

(143) 


We point out that the Coulomb interaction is not well 
suited for the fc-constant coupling approximation. The 
reason is that the function 9 f l„(J) behaves as 1/n for large 
n while Z n {J) goes to 1 and thus the Matsubara integral 
shows a logarithmic divergence due to C atat (0, J, 0, J') 
if A^(J) does not cut off the integral. Often the ef¬ 
fect of the Coulomb potential is mimicked by replacing 
C atat with /x*0(w c - |w n |) where = 1+Cst S & i ° (£/uJc ) with 
£, a parameter of the electronic band structure and w c 
a phonon frequency cutoff 61,62 . Usually the so called 
Morel-Anderson pseudo potential /r * is fitted so that the 
calculated T c matches the experimental one. fi* usually 
ranges between 0.1 and 0.16 for conventional SC 2 '. The 
above equations imply that the coupling is isotropic in 
the sense that all states with equal center of energy and 
equal splitting share the same coupling matrix elements. 
Sometimes as in the well known case of MgB 2 there are 
significant differences in the couplings and it is important 
to group states into bands for the isotropic approximation 
to hold. We refer to this case as the multiband approx¬ 
imation which simply means that all isotropic variables 
obtain another index for the band they correspond to. 

Comparing the equations for the SC KS GF of Eq. (65) 
(noting u££(«k*")* = asign(cr)A a fc /F fe where F k = 

^ e kr +e. ki + |^ fc | 2 ) with the interacting GF Eq. (123) 

we note that Z\ a fe takes the role of A k (oj n ) so the similar 
name is not accidental. However, as we have seen A k (u} n ) 
takes its significant shape in Matsubara space while A s sk 
does not have such a dependence and mimics the SC 
pairing in its k dependence in a way that densities of the 
interacting system are reproduced. 


V. SUMMARY AND CONCLUSION 

In this work we have developed fully ab-initio methods 
to compute the SC phase of a material in a magnetic field 
Zeeman-coupled to the spin magnetization. In a unified 
notation we present a purely GF based (the Eliashberg 
approach) and a Density Functional based scheme. 

In our DFT we have employed a SC KS system to 
reproduce the interacting densities tt(t’), m(r), x(r, r') 
and r(Ri... Rn). The SC KS system can be solved 
analytically using the SDA where we only consider the 
singlet pairing of time reversed basis states. We have de¬ 
rived £C-potentials in this case that include the electron- 
nuclear interaction on the level of KS phonons and treats 
the Coulomb interaction in the same footing without the 
need for any adjustable parameter. 
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As a second step we have applied similar approxima¬ 
tions to the Dyson equation starting from the SC KS 
system as a formally non-interacting system. This pro¬ 
cedure leads to the Eliashberg equations of a SC in a 
magnetic field similar to those discussed in Ref. 24. 

While SpinSCDFT allows to include the full Coulomb 
potential and promises numerically efficient calculations 
tor real materials, the direct GF approach is, instead, 
valuable to get direct physical insights to develop approx¬ 
imations and further improve the SpinSCDFT scheme. 

The theoretical framework presented in this work al¬ 
lows to compute the phenomenon of coexistence and com¬ 
petition of SC with magnetism from first principles. Es¬ 
pecially in connection with the discovery of Fe supercon¬ 
ductors this was intensively studied in recent years. 

In the subsequent part II, we will discuss a detailed 
numerical implementation of the equations presented in 
this work, i.e. the linear and non-linear functionals and 
the Eliashberg equations without Coulomb interactions. 
Further we will introduce a GoWo scheme to obtain the 
excitation spectrum starting from a SpinSCDFT calcu¬ 
lation. 


Appendix A: Formulas For The Matsubara Sums 


In the potential terms it appears the Matsubara sum¬ 
mation 


r,(£,E') = s Es— 


P „ On - E)(iu n - E') 


(Al) 


This is analytically evaluated with the result 


Ps(E,E') 


fp{E) — fp(E') 

E~E’ 


(A2) 


lim P S (E,E’) = d E fp(E) = -fifp(E)fp(-E)(A3) 

h/ — yh/ 


where the symmetries P S (E,E') = P s (—E,—E') and 
P S (E,E') = P S (E',E) hold. The Matsubara frequency 
summation 


I(fl 1 E 1 ,E 2 ,E 3 )-—^2- 


P 2 O ~ O i K - - ft 


1 


1 


p p ( A4 ) 

\U> n ' — A 2 lU> n — Jh 3 

L(ft, E u E 2 , E 3 ) = I(-ft , E u E 2 , E 3 ) -I(ft, E u E 2 , E 3 ) 

(A5) 


is also in principle straightforward. However the result¬ 
ing formulas are rather large and computer algebra be¬ 
comes essential for the evaluation of residues and limiting 
behaviours, necessary for a numerical implementation. 
Note that a partial summation leads to 


1 


L(ft,E 1 ,E 2 ,E 3 ) = -J2 


M pb (ft, E 2 , U>n) 


, , w ,. (A6) 

P ^ (lU n - £i) (lW„ - E 3 ) 

From the definition we observe the following symmetry 
relations 


L(ft,E 1 ,E 2 ,E 3 ) = L{ft,E 3 ,E 2 ,E 1 ) (A7) 

L{-ft, Ex, E 2 , E 3 ) = -L(ft, E 3 , E 2 , Ex) (A8) 
(L(ft, E u E 2 , E 3 ))* = -L(ft, -E x , -E 2 , -E 3 ) (A9) 

Evaluation of the Coulomb requires the following sum¬ 
mation 

L a (E 1 ,E 2 ,E 3 ) = ±Y liu}n ,\ E2 x 

nn' 

X . 1 p ■ 1 p (A10) 

lU>n — -Cu 1 iO n — ht 3 

= fp{E 2 )P B (E 1 ,E 3 ) (All) 

Using Mathematica, we evaluate the sums Eqs. (A4) and 
(A5) to 


L(HE E E)-( f^Mft) | fp(E 2 )(l + np(ft)) f p (Ex){l-fp(E 2 ) + np(C2)) 

{ ’ U 2 ’ 31 y (E 2 — Ei + tl)(E 2 -E 3 + n) (Ex -E 2 + ft)(E 3 -E 2 + ft) (Ex - E 3 )(Ex — E 2 — ft) 

fp(E 3 )(l-f 0 (E 2 )+np(n)) f p (Ex)(f 0 (E 2 )+np(E!)) fp(E 3 ){f p (E 2 ) +n p (ft)) \ 

(Ex - E 3 )(E 2 - E 3 + Q) + (Ex-E 3 )(Ex-E 2 + ft) + (E 3 -Ex)(E 3 -E 2 + ft)J' [ j 


Clearly some points, e.g. Ex = E 3 are numerically problematic, so whenever Ex ~ E 3 we may have to evaluate the 
limiting formula instead. In general, the various limits where the denominators are zero, all exist and can be computed 
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explicitly, again using Mathematica. The results are 


lim 

E\—¥E% 


L(f2,E 1 ,E 2 ,E 3 ) = fp(E 2 ) 


fp{E 2 ) + np(fi) _ 

(e 2 -e 3 - ny {e 2 


( np{f2) 1 +rif}(n) \ _ / fp(-E 2 ) + np(Q) 

\{E 2 ~E 3 + ny ( E 2 -E 3 -ny) m 3J V {E 2 -E 3 + ny 

^ f E 3 y E -\ny + ( 1 2 3 4 5 6 * 8 9 * 11 12 Mfl) +i)(^ - ^))) , (Ai3) 


lim lim L(f2,Ei,E 2 ,E 3 ) = B 

Q^E 3 -E 2 -Ei->_E 3 


(1 + fp(E 2 ) + np{E 3 - E 2 ))fp(-E 3 )fp(E 3 ) 


fp(E 2 ) + fp(E 3 )(l~2fp(E 2 )) 
4{E 2 -E 3 y 


2{E 2 - E 3 ) 

+ j3 2 fp(—E 2 )(2 + np(E 3 - E 2 ))fp(E 3 )(± - fp(E 3 )) , (A14) 


_ lim „ L(f2,E 1 ,E 2 ,E 3 ) 

1 1—y E i —E 2 


/ /3 (£ 1 )(/ /3 (£ 2 )+n /3 (£; 1 -.E; 2 )) / /3 (^2)(l + ^(£ 1 -£ 2 )) 

2 (E! - E 2 ){E 1 - E 3 ) + 2(£7 1 - ^ 2 )( j B 1 - 2E 2 + E 3 ) + 

fp{E 3 )(np(E 1 - E 2 ) + fp{-E 2 )) - fp{E 2 )np{E 1 - E 2 ) 

(E ± - E 3 y 

fp{E 3 )(fp{E 2 ) + np{E 1 - E 2 )) M-EJfpj^npfa-Ei) 
2(E 3 -E 1 )(E 1 -2E 2 + E 3 ) P E 3 -E! 


(A15) 


lim lim L(f2 1 E\, E 2 , E 3 ) = 

E 1 ^2E 2 -E 3 q^e 1 -e 2 


fp{E 2 )(l + np{E 2 ~E 3 )) 
2 (E 2 — E 3 ) 


(3fp(-E 3 ) - 


1 


2 {E 2 - E 3 ) 


(A16) 


We point out here that the Limit Q q \ —>■ 0 does not exist. It is however unimportant as the go to zero in the 
limit fl —> 0 faster than L diverges. 
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